

program cs307_matrix_calculator
implicit none

integer, parameter :: NMAX = 500
integer::fd0=0,fd1=1,fd2=2,fd3=3,st0,st1,st2,NOM,result_nnz,info,s,errorflag,m,n,in_row,in_col
integer::nnz_1,rowNum_1,ColNum_1,i,j,k,p,nnz_2,rowNum_2,Colnum_2,it
double precision::result,DMGT,eps,det,d1,d2
character(50)::CommandFile,MatrixFile1,MatrixFile2
character(1)::OType
character(10)::rep
character(19)::sym
character(7)::field

integer,allocatable,dimension(:):: Ai_1, Aj_1, ia_1, ja_1,Ai_2,Aj_2,ia_2,ja_2,result_ia,result_ja
double precision,allocatable,dimension(:):: Aa_1,a_1,Aa_2,a_2,dummy,result_a,temp,vector_result
double precision,allocatable,dimension(:,:)::matrix_1,matrix_2,matrix_result,dummy2,A


!set all the parameter and flags
CommandFile="command.txt"

!read in the file
open (fd0,file=CommandFile,status='old',action='read',iostat=st0)

if(st0 .ne. 0) then
   print *, "the command file read fails, program terminated"
   stop
endif

!read in the operation type
read(fd0,"(A1)") OType


!print *,OType


!check the operation type to determine if it is valid and determine the number of the matrices
if( (OType .eq. "A") .or. (OType .eq. "S") .or. (OType .eq. "M") .or. (OType .eq. "V") .or. (OType .eq. "O") .or. (OType .eq. "Z"))then
   NOM = 2
!   print *, "the NOM is set to 2"
else if ((OType .eq. "D") .or. (OType .eq. "T") .or. (OType .eq. "I") .or. (OType .eq. "E") .or. (OType .eq. "C")) then
   NOM = 1
!   print *, "the NOM is set to 1"
else 
   print *,"The operation type is invalid"
   stop
endif

!open result.txt, which will store the information about the matrix
open (fd3,file="result.mtx",action='write',iostat=st1)

if(st1 .ne. 0) then
   print *, " result.mtx cannot be opened, program terminated"
   stop
endif




!read in the path of the taget matrix/matrices
read(fd0,*) MatrixFile1
!   print *, MatrixFile1

!open matrix 1
open (fd1,file=MatrixFile1,status='old',action='read',iostat=st1)

if(st1 .ne. 0) then
   print *, " matrix1 read fails, program terminated"
   stop
endif




!get info and read in matrix if requested
if(NOM .eq. 2) then
   read(fd0,*) MatrixFile2
!   print *, MatrixFile2
   open (fd2,file=MatrixFile2,status='old',action='read',iostat=st2)

   if(st2 .ne. 0) then
      print *, " matrix2 read fails, program terminated"
      stop
   endif

endif





!read in the actually matrix1
call mminfo(fd1,rep,field,sym,rowNum_1,colNum_1,nnz_1)

allocate( Aa_1( nnz_1), Ai_1( nnz_1), Aj_1( nnz_1),dummy( nnz_1))
call mmread(fd1 ,rep,field,sym, rowNum_1, colNum_1, nnz_1, nnz_1, Ai_1, Aj_1,dummy, Aa_1,dummy)
deallocate(dummy)

!map it to a real matrix
allocate(matrix_1(rowNum_1,colNum_1))
matrix_1 = 0
do i=1,nnz_1
   matrix_1(Ai_1(i),Aj_1(i))=Aa_1(i)
enddo



allocate( a_1( nnz_1), ia_1( rowNum_1+1 ), ja_1( nnz_1))
call coocsr( rowNum_1, nnz_1, Aa_1, Ai_1, Aj_1, a_1, ja_1, ia_1)
deallocate( Aa_1, Ai_1, Aj_1)







!print *, a_1,ja_1,ia_1

!read in matrix 2 if applicable
if(NOM .eq. 2)then

   
   call mminfo(fd2,rep,field,sym,rowNum_2,colNum_2,nnz_2)

   allocate( Aa_2( nnz_2), Ai_2( nnz_2), Aj_2( nnz_2),dummy( nnz_2))
   call mmread(fd2 ,rep,field,sym, rowNum_2, colNum_2, nnz_2, nnz_2, Ai_2, Aj_2,dummy, Aa_2,dummy)
   deallocate(dummy)

   allocate(matrix_2(rowNum_2,colNum_2))
   matrix_2 = 0
   do i=1,nnz_2
      matrix_2(Ai_2(i),Aj_2(i))=Aa_2(i)
   enddo


   allocate( a_2( nnz_2), ia_2( rowNum_2+1 ), ja_2( nnz_2))
   call coocsr( rowNum_2, nnz_2, Aa_2, Ai_2, Aj_2, a_2, ja_2, ia_2)
   deallocate( Aa_2, Ai_2, Aj_2)


endif







!addition operation
if( OType .eq. "A"  ) then

   if( (rowNum_1 .ne. rowNum_2) .or. (colNum_1 .ne. colNum_2)) then
      print *,"the two matrices dimension do not agree, cannot perform addition, program terminated"
      stop
   endif

   allocate(matrix_result(rowNum_1,colNum_1))
   matrix_result = 0
   !$OMP PARALLEL DO PRIVATE(i,j)
   do i=1,rowNum_1
      do j=1,ColNum_1
         matrix_result(i,j) = Matrix_1(i,j) + Matrix_2(i,j)
      enddo
   enddo
 !$OMP END PARALLEL DO
   write(fd3,*)'%%MatrixMarket matrix coordinate real general'
   write(fd3,*)rowNum_1,colNum_1,colNum_1*rowNum_1

   do i=1,rowNum_1
      do j=1,ColNum_1
         write(fd3,*) i,j,matrix_result(i,j)
      enddo
   enddo

   close(fd3)
   print *, 'READY'
   stop
   

endif






!substraction operation
if(OType .eq. "S" )then

   if( (rowNum_1 .ne. rowNum_2) .or. (colNum_1 .ne. colNum_2)) then
      print *,"the two matrices dimension do not agree, cannot perform subtraction, program terminated"
      stop
   endif

   allocate(matrix_result(rowNum_1,colNum_1))
   matrix_result = 0
!$OMP PARALLEL DO PRIVATE(i,j)
   do i=1,rowNum_1
      do j=1,ColNum_1
         matrix_result(i,j) = Matrix_1(i,j) - Matrix_2(i,j)
      enddo
   enddo
 !$OMP END PARALLEL DO
   write(fd3,*)'%%MatrixMarket matrix coordinate real general'
   write(fd3,*)rowNum_1,colNum_1,colNum_1*rowNum_1

   do i=1,rowNum_1
      do j=1,ColNum_1
         write(fd3,*) i,j,matrix_result(i,j)
      enddo
   enddo

   close(fd3)
   print *, 'READY'
   stop


endif




!mulitplication with another matrix or vector
if(OType .eq. "M" .or. OType .eq. "V")then
   
   if(  rowNum_2 .ne. colNum_1) then
      print *,"the two matrices dimension do not agree, cannot perform multiplication, program terminated"
      stop
   endif

   allocate(matrix_result(rowNum_1,colNum_2))
   matrix_result = 0
   !$OMP PARALLEL DO PRIVATE(s,i,k,j)
   do i=1,rowNum_1
      do k=ia_1(i),ia_1(i+1)-1
         j=ja_1(k)
         do s=1,colNum_2
            matrix_result(i,s)=matrix_result(i,s)+a_1(k)*matrix_2(j,s)
         enddo
      enddo
   enddo
   !$OMP END PARALLEL DO
   write(fd3,*)'%%MatrixMarket matrix coordinate real general'
   write(fd3,*)rowNum_1,colNum_2,colNum_1*rowNum_1

   do i=1,rowNum_1
      do j=1,ColNum_2
         write(fd3,*) i,j,matrix_result(i,j)
      enddo
   enddo

   close(fd3)
   print *, 'READY'
   stop
   
   
   !!write the result matrix to file(result.mtx)
   



endif



!determinant calculation
if(OType .eq. "D") then

   if(rowNum_1 .ne. colNum_1) then
      print *, "the matrix is not square, cannot performe determinant calculation! Program terminated!"
      stop
   endif

   result= DMGT(eps,rowNum_1,matrix_1)
   
   write(fd3,*)'%%MatrixMarket matrix coordinate real general'
   write(fd3,*)"1","1","1"

   write(fd3,*)"1","1",result

   close(fd3)
   print *, 'READY'
   stop
   
 
   
endif


!transpose calculation
if(OType .eq. "T") then
   
   allocate(matrix_result(colNum_1,rowNum_1))
   !$OMP PARALLEL DO PRIVATE(i,j)
   do i=1,rowNum_1
      do j=1,colNum_1
         matrix_result(i,j) = matrix_1(j,i)
      enddo
   enddo
 !$OMP END PARALLEL DO
   write(fd3,*)'%%MatrixMarket matrix coordinate real general'
   write(fd3,*)colNum_1,rowNum_1,colNum_1*rowNum_1
   do j=1,ColNum_1
      do i=1,rowNum_1
     	write(fd3,*) i,j,matrix_result(i,j)
      enddo
   enddo

   close(fd3)
   print *, 'READY'
   stop
   
   
   print *, matrix_result

endif




!!this solver can only solve for A*X=B where A is n*n and B is n*m
if(OType .eq. "O") then
  
   if((rowNum_1 .ne. colNum_1) .or. (colNum_1 .ne. rowNum_2)) then
      print *, "the dimension do not agree to do a solving, program terminated!"
      stop
   endif
   
   allocate(matrix_result(colNum_1,colNum_2))
   
   call MATINV(rowNum_1,colNum_2,matrix_1,matrix_2,det)
   
   write(fd3,*)'%%MatrixMarket matrix coordinate real general'
   write(fd3,*)colNum_1,colNum_2,colNum_1*colNum_2

   do i=1,colNum_1
      do j=1,colNum_2
         write(fd3,*) i,j,matrix_result(i,j)
      enddo
   enddo

   close(fd3)
   print *, 'READY'
   stop
   
endif


!inverse calculation
if(OType .eq. "I") then
   
   if(rowNum_1 .ne. colNum_1) then
      print *, "the matrix is not square, cannot performe inverse calculation! Program terminated!"
      stop
   endif

   allocate(matrix_result(rowNum_1,colNum_1))
   allocate(dummy2(rowNum_1,1))

   !!mkl routine
   dummy2=0
   do i=1,rowNum_1
      do j=1,colNum_1
         dummy2(i,1) = dummy2(i,1)+matrix_1(i,j)
      enddo
   enddo
   
   call MATINV(rowNum_1,1,matrix_1,dummy2,det)
    
   write(fd3,*)'%%MatrixMarket matrix coordinate real general'
   write(fd3,*)rowNum_1,colNum_1,colNum_1*rowNum_1

   do i=1,rowNum_1
      do j=1,ColNum_1
         write(fd3,*) i,j,matrix_result(i,j)
      enddo
   enddo

   close(fd3)
   print *, 'READY'
   stop
   
   
   
   
   
endif





!compute the eigenvalue of a matrix, matrix has to be smaller than 500*500
if(OType .eq. "E") then
   if(rowNum_1 .ne. colNum_1) then
      print *, "the matrix is not square, cannot performe eigenvalue computing calculation! Program terminated!"
      stop
   endif
   
   if(rowNum_1 .gt. 500) then
      print *,"target matrix is too big to compute eigenvalue"
      stop
   endif
   
   eps=1.d-10
   d1=1.d-8
   d2=1.d-8
   m = 600
   n = rowNum_1
   allocate(vector_result(rowNum_1))
   allocate(matrix_result(rowNum_1,rowNum_1))
   
   allocate(A(NMAX,NMAX))
   A = 0
   do i=1,rowNum_1
      do j=1,colNum_1
         A(i,j) = matrix_1(i,j)
      enddo
   enddo

      
   call EPMRI(eps,d1,d2,m,n,A,IT,vector_result,matrix_result)

   !the eigenvalue is stored in the vector_result
   if(IT .eq. -1) then
      print *, "warning! no convergence is achieved, the result might not be accurate!"
      stop
   endif
 

   write(fd3,*)'%%MatrixMarket matrix coordinate real general'
   write(fd3,*)rowNum_1,"1",rowNum_1

   do i=1,rowNum_1
      
         write(fd3,*) i,"1",vector_result(i)
     
   enddo

   close(fd3)
   print *, 'READY'
   stop
   
   
   

endif





!compute the eigenvector of a matrix, matrix has to be smaller than 500*500
if(OType .eq. "C") then

   if(rowNum_1 .ne. colNum_1) then
      print *, "the matrix is not square, cannot performe eigenvector computing calculation! Program terminated!"
      stop
   endif


   if(rowNum_1 .gt. 500) then
      print *,"target matrix is too big to compute eigenvector"
      stop
   endif

   eps = 1.d-10
   d1 = 1.d-8
   d2 = 1.d-8
   m = 600
   n = rowNum_1

   allocate(vector_result(rowNum_1))
   allocate(matrix_result(rowNum_1,rowNum_1))

   vector_result = 0
   matrix_result = 0

   allocate(A(NMAX,NMAX))
   A = 0
   do i=1,rowNum_1
      do j=1,colNum_1
         A(i,j) = matrix_1(i,j)
      enddo
   enddo

   call EPMRI(eps,d1,d2,m,n,matrix_1,IT,vector_result,matrix_result)

   if(IT .eq. -1) then
      print *, "warning! no convergence is achieved, the result might not be accurate!"
      stop
   endif
   !the eigenvector is stored in the matrix_result

       write(fd3,*)'%%MatrixMarket matrix coordinate real general'
   write(fd3,*)rowNum_1,rowNum_1,rowNum_1*rowNum_1

   do i=1,rowNum_1
      do j=1,rowNum_1
         write(fd3,*) i,j,matrix_result(i,j)
      enddo
   enddo

   close(fd3)
   print *, 'READY'
   stop
   
   
endif


!The user like to define their own procedure
!if(OType .eq. "Z") then


!since the operation type is self-defined, there is no dimension check for the matrices.
!the interface is followed the format Subroutine SelfDefined(matrix_1,matrix_2,rowNum_1,colNum_1,rowNum_2,colNum_2,result,matrix_result)
!   read(fd0,*) in_row
!   read(fd0,*) in_col
   
!   allocate(matrix_result(in_row,in_col))
   
!   call SelfDefined(matrix_1,matrix_2,rowNum_1,colNum_1,rowNum_2,colNum_2,matrix_result)
   
!   print *,"according to the user defined subroutine, the matrix_result is ",matrix_result

!endif



end program


!The function DMGT returns the determinant of a real square matrix
!A(n,n) by Gauss method with full pivoting.
!----------------------------------------------------------------------------
!  Input parameters:
!  eps        precision (real*8)
!  n          size of A matrix (integer)
!  A          pointer to input real square matrix
!  Output parameters:
!  None
!-----------------------------------------------------------------------------
!The procedure TSRGT is used to reduce A matrix to an upper triangular matrix.
!Output variables are it(integer), C(n,n), Kp(n) and Lp(n).
!If it=0, matrix A is singular, if it=1, matrix A is regular. Table C contains
!at location i,j (j>=i) the corresponding element of the upper triangular matrix.
!Tables Lp and Kp contain informations relative to exchanges of line or column
!that occured during the process. For instance, the element number k of Lp is
!an integer <> k if an exchange of line has been made at step k (k=1,2,...,n).
!The number of exchanges of lines and columns is stored in integer L. the
!determinant of A matrix is stored in d0 (real*8).
!-----------------------------------------------------------------------------
real*8 Function DMGT(eps, n, A)
integer n
real*8 eps, A(n,n)
real*8 d0

real*8, pointer :: C(:,:)
integer,pointer :: Kp(:), Lp(:)

!allocate local matrix C and vectors Kp, Lp
  allocate(C(n,n),STAT=ialloc)
  allocate(Kp(n),STAT=ialloc)
  allocate(Lp(n),STAT=ialloc)

  call TSRGT(eps,n,A,it,C,Kp,Lp)  !call triangularization subroutine
  if (it==0) then
    d0=0.d0  !matrix singular, det=0
  else       !matrix regular, det<>0
    d0=1.d0
    do k=1, n
	  d0=d0*C(k,k)
    end do
    l=0
    do k=1, n-1
      if (Lp(k).ne.k)  l=l+1
      if (Kp(k).ne.k)  l=l+1
    end do
    if (MOD(l,2).ne.0) d0=-d0  !l is odd
  end if
  DMGT=d0   !return determinant
  return
End

!End of file deter.f90


!The subroutine TSRGT applies to input real square matrix A(n,n) the upper
!triangularization algorithm of Gauss method with full pivoting and keeps
!trace of successive transformations done in integer vectors KP and LP.
!-----------------------------------------------------------------------------
!  Input parameters:
!  eps        precision (real*8)
!  n          size of A matrix (integer)
!  A          pointer to input real square matrix (real*8)
!  Output parameters:
!  it         flag=1 if A matrix ok, =0 if A matrix is singular (integer)
!  C          pointer to table storing main diagonal elements and supra-
!             diagonal elements of upper triangular matrix and the multi-
!             plying coefficients used during triangularization process
!  KP         table storing informations concerning the column exchanges
!             during process (integer)
!  LP         table storing informations concerning the line exchanges
!             during process (integer)
!-----------------------------------------------------------------------------
!The table C is first initialized to A matrix, then receives at each step k
!of the triangularization process, usefull elements of A matrix at step k for
!k=1,2,...n.
!The variables po(real*8), lo and ko(integer) store respectively pivot at step k,
!its line number and its column number.
!------------------------------------------------------------------------------
Subroutine TSRGT(eps, n, A, it, C, Kp, Lp)
  real*8 eps
  integer n,it
  real*8 A(n,n), C(n,n)
  integer Kp(n),Lp(n) 
  real*8  po,t0
  C=A; it=1; k=1
  do while (it==1.and.k<n)
    po=C(k,k); lo=k; ko=k
    do i=k, n
      do j=k, n
        if (dabs(C(i,j))>dabs(po)) then
          po=C(i,j); lo=i; ko=j
        end if
      end do
    end do
    Lp(k)=lo; Kp(k)=ko
    if (dabs(po)<eps) then
      it=0
    else
      if (lo.ne.k) then
        do j=k, n
          t0=C(k,j); C(k,j)=C(lo,j); C(lo,j)=t0
        end do
      end if
      if (ko.ne.k) then
        do i=1, n
          t0=C(i,k); C(i,k)=C(i,ko); C(i,ko)=t0
        end do
      end if 
      do i=k+1, n
        C(i,k)=C(i,k)/po
        do j=k+1, n
          C(i,j)=C(i,j)-C(i,k)*C(k,j)
        end do 
      end do
      k=k+1
    end if
  end do
  if (it==1.and.dabs(C(n,n))<eps)  it=0
  return
End !TSRGT






!*******************************************
!*  SOLVING A LINEAR MATRIX SYSTEM AX = B  *
!*  with Gauss-Jordan method using full    *
!*  pivoting at each step. During the pro- *
!* cess, original A and B matrices are     *
!* destroyed to spare storage location.    *
!* --------------------------------------- *
!* INPUTS:    A   MATRIX N*N               *
!*            B   MATRIX N*M               *
!* --------------------------------------- *
!* OUTPUTS:   A   INVERSE OF A N*N         *
!*            DET  DETERMINANT OF A        *
!*            B   SOLUTION MATRIX N*M      *
!* --------------------------------------- *
!* NOTA - If M=0 inversion of A matrix     *
!*        only.                            *
!*******************************************
      SUBROUTINE MATINV(N,M,AA,BB,DET)                                
      REAL*8, PARAMETER :: EPSMACH=2.D-16                                    
      REAL*8 AA(N,N),BB(N,M)                                                   
      REAL*8,POINTER ::  PC(:),PL(:),CS(:)                                                  
      REAL*8 PV,PAV,DET,TT    
	                                                    
!Initializations :                       
      allocate(PC(1:N),stat=ialloc)
      allocate(PL(1:N),stat=ialloc)
      allocate(CS(1:N),stat=ialloc)	                                        
      DET=1.D0
      DO I=1,N                                                                                                                   
        PC(I)=0.D0
	PL(I)=0.D0
	CS(I)=0.D0
      END DO
!main loop                                                                   
      DO K=1,N                                                                  
!Searching greatest pivot :                                               
        PV=AA(K,K)                                                              
        IK=K                                                                    
        JK=K                                                                    
        PAV=DABS(PV)                                                            
        DO I=K,N                                                                
          DO J=K,N                                                              
            IF (DABS(AA(I,J)).GT.PAV) THEN                                      
              PV=AA(I,J)                                                        
              PAV=DABS(PV)                                                      
              IK=I                                                              
              JK=J                                                              
            ENDIF                                                               
          ENDDO                                                                 
        ENDDO                                                                   
                                                                               
!Search terminated, the pivot is in location I=IK, J=JK.
!Memorizing pivot location: :                                        
        PC(K)=JK                                                                
        PL(K)=IK                                                                
                                                                               
!Determinant DET is actualised
!If DET=0, ERROR MESSAGE and STOP
!Machine dependant EPSMACH equals here 1.DE-20                                        
                                                                               
        IF (IK.NE.K) DET=-DET                                                   
        IF (JK.NE.K) DET=-DET                                                   
        DET=DET*PV                                                              
        IF (DABS(DET).LT.EPSMACH) THEN                                          
!Error message and Stop                                                     
          PRINT 10                                                              
          STOP                                                                  
        ENDIF                                                                   
                                                                               
!POSITIONNING PIVOT IN K,K:                                              
        IF(IK.NE.K) THEN                                                        
          DO I=1,N                                                              
!EXCHANGE LINES IK and K:                                            
            TT=AA(IK,I)                                                         
            AA(IK,I)=AA(K,I)                                                    
            AA(K,I)=TT                                                          
          ENDDO                                                                 
        ENDIF                                                                   
      IF (M.NE.0) THEN                                                          
        DO I=1,M                                                               
          TT=BB(IK,I)                                                           
          BB(IK,I)=BB(K,I)                                                      
          BB(K,I)=TT                                                            
        ENDDO                                                                   
      ENDIF                                                                     
!Pivot is at correct line                                                
        IF(JK.NE.K) THEN                                                        
          DO I=1,N                                                              
!Exchange columns JK and K of matrix AA                                         
            TT=AA(I,JK)                                                         
            AA(I,JK)=AA(I,K)                                                    
            AA(I,K)=TT                                                          
          ENDDO                                                                 
        ENDIF                                                                   
!Pivot is at correct column and located in K,K                                              
                                                                               
!Store column K in vector CS                             
!then set column K to zero                                             
        DO I=1,N                                                                
          CS(I)=AA(I,K)                                                         
          AA(I,K)=0.D0                                                          
        ENDDO                                                                   
!                                                                               
        CS(K)=0.                                                                
        AA(K,K)=1.                                                              
!Modify line K :                                            
        IF(DABS(PV).LT.EPSMACH) THEN                                            
          WRITE(*,*) '  PIVOT TOO SMALL - STOP'                               
          STOP                                                                  
        ENDIF                                                                   
        DO I=1,N                                                                
          AA(K,I)=AA(K,I)/PV                                                    
        ENDDO                                                                   
        IF (M.NE.0) THEN                                                        
          DO I=1,M                                                             
            BB(K,I)=BB(K,I)/PV                                                  
          ENDDO                                                                 
        ENDIF                                                                   
!Modify other lines of matrix AA:                                        
        DO J=1,N                                                                
          IF (J.EQ.K) CONTINUE                                                  
          DO I=1,N                                                              
!Modify line J of matrix AA :                                            
            AA(J,I)=AA(J,I)-CS(J)*AA(K,I)                                       
          ENDDO                                                                 
          IF (M.NE.0) THEN                                                      
            DO I=1,M                                                          
              BB(J,I)=BB(J,I)-CS(J)*BB(K,I)                                     
            ENDDO                                                               
          ENDIF                                                                 
        ENDDO                                                                   
!Line K is ready.                                                
      ENDDO                                                                     
!End of K loop                                                              
                                                                               
!The matrix AA is inverted - Rearrange AA                         
                                                                               
!Exchange lines                                                            
      DO I=N,1,-1                                                               
        IK=PC(I)                                                                
        IF (IK.EQ.I) CONTINUE                                                   
!EXCHANGE LINES I AND PC(I) OF AA:                                         
        DO J=1,N                                                                
          TT=AA(I,J)                                                            
          AA(I,J)=AA(IK,J)                                                      
          AA(IK,J)=TT                                                           
        ENDDO                                                                   
        IF (M.NE.0) THEN                                                        
          DO J=1,M                                                             
            TT=BB(I,J)                                                          
            BB(I,J)=BB(IK,J)                                                    
            BB(IK,J)=TT                                                         
          ENDDO                                                                 
        ENDIF                                                                   
!NO MORE EXCHANGE NEEDED                                                      
!GO TO NEXT LINE                                                  
      ENDDO                                                                     
                                                                               
!EXCHANGE COLUMNS                                                          
                                                                               
      DO J=N,1,-1                                                               
        JK=PL(J)                                                                
        IF (JK.EQ.J) CONTINUE                                                   
!EXCHANGE COLUMNS J AND PL(J) OF AA :                                       
        DO I=1,N                                                                
          TT=AA(I,J)                                                            
          AA(I,J)=AA(I,JK)                                                      
          AA(I,JK)=TT                                                           
        ENDDO                                                                   
!NO MORE EXCHANGE NEEDED                                                      
!GO TO NEXT COLUMN   
      ENDDO                                                                     
!REARRANGEMENT TERMINATED.                                                        
      RETURN                                                                    
   10 FORMAT(///'  DETERMINANT EQUALS ZERO, NO SOLUTION!')                    
      END          



                                                             

  !***************************************************************
  !* Subroutine DECCRM determines the lower triangular matrix and*
  !* the upper triangukar matrix of Crout's decomposition of a   *
  !* given square real matrix, A.                                *
  !* ----------------------------------------------------------- *
  !* INPUTS:                                                     *
  !*         eps: required precision (double)                    *
  !*          n : size of matrix A (integer)                     *
  !*          A : input matrix (n x n)                           *
  !* OUTPUTS:                                                    *
  !*          it: flag, =0 if method does not apply              *
  !*                    =1 if method is ok.                      *
  !*           U: lower triangular matrix.                       *
  !*           V: upper triangular matrix.                       *
  !***************************************************************
  Subroutine DECCRM(eps, n, A, it, U, V)
    integer, parameter :: NMAX = 500
    real*8 A(NMAX,NMAX), U(NMAX,NMAX), V(NMAX,NMAX)
    real*8 eps, s


   
    if (dabs(A(1,1)) < eps)  then
	  it=0
    else
      do i=1, n
	    U(i,1)=A(i,1)
      end do 
      V(1,1)=1.d0
      do j=2, n
	    V(1,j)=A(1,j)/U(1,1)
      end do
      it=1; k=2
      do while (it.ne.0.and.k<=n)
	    do i=1, n
	      if (i < k) then  
		    U(i,k)=0.d0
		  else
	        s=0.d0
	        do j=1, k-1
			  s = s + U(i,j)*V(j,k)
            end do
	        U(i,k) = A(i,k) - s
		  end if
        end do
	    if (dabs(U(k,k)) < eps) then 
          it=0
		else
	      do j=1, n
	        if (j < k) then 
			  V(k,j)=0.d0
	        else if (j==k) then
			  V(k,j)=1.d0
            else
	          s=0.d0
	          do i=1, k-1
			    s = s + U(k,i)*V(i,j)
              end do
	          V(k,j) = A(k,j)/U(k,k);
			end if
          end do
	      k = k + 1
		end if
      end do !while
    end if
	return
  End

!  void MatPrint(char *title, MAT A, int n) {
!  int i,j;
!   printf("\n %s\n", title);
!for (i=1; i<=n; i++) {
!     for (j=1; j<=n; j++)
!       printf(" %f", A[i][j]);
!	  printf("\n");
!   }
! }

! void VecPrint(char *title, VEC X, int n) {
!   int i;
!   printf("\n %s\n", title);
!   for (i=1; i<=n; i++) printf(" %f", X[i]);
!	printf("\n");
! }

  !*********************************************************
  !* Calculate the eigenvalues of a real square matrix by  *
  !* Rutishauser's Method.                                 *
  !* ----------------------------------------------------- *
  !* INPUTS:                                               *
  !*        eps: absolute precision (double)               *
  !*        dta: relative precision (double)               *
  !*         m : maximum number of iterations (integer)    *
  !*         n : size of matrix A (integer)                *
  !*         A : input real square matrix (n x n)          *
  !* OUTPUTS:                                              *
  !*         it: flag, =-1 if convergence is not obtained  *
  !*                   =1 if convergence is ok.            *
  !*         R : contains in output the n eigenvalues of A *
  !*                                                       *         
  !*********************************************************
  Subroutine VAMR(eps, dta, m, n, A, it, R)
    integer, parameter :: NMAX = 500
    real*8 A(NMAX,NMAX), R(NMAX)
    real*8 eps, dta
    real*8 phi,s,t0
    real*8 U(NMAX,NMAX), V(NMAX,NMAX)

    t0=0.d0
    l=1
    do while (l<=m.and.it.ne.1)
      do i=1, n
	    R(i)=A(i,i)
      end do

      call DECCRM(eps, n, A, it, U, V)

      if (it==0) then
	    do i=1, n
		  A(i,i)=A(i,i) + 1.d0
        end do
	    t0 = t0 + 1.d0
      else
	    do i=1, n
		  do j=1, n
	        s=0.d0
	        do k=1, n
			  s = s + V(i,k) * U(k,j)
            end do
	        A(i,j) = s
		  end do
        end do
	    phi=0.d0
        do i=1, n
	      s= dabs(A(i,i)-R(i))
	      if (s > phi)  phi=s
		end do
	    if (phi < dta) then
	      do i=1, n
		    R(i) = A(i,i) - t0
          end do
		else
	      l = l + 1
	      it=-1
		end if
      end if
    end do !while
	return
  End

  !************************************************************
  !* Procedure IIM calculates a real eigenvalue and the asso- *
  !* ciated eigenvector of a real square matrix the inverse   *
  !* iteration method.                                        *
  !* -------------------------------------------------------- *
  !* INPUTS:                                                  *
  !*         eps : absolute precision (double)                *
  !*         dta : relative precision (double)                *
  !*          m  : maximum number of iterations (integer)     *
  !*          n  : size of matrix A                           *
  !*          A  : input real square matrix (n x n)           *
  !* OUTPUTS:                                                 *
  !*          it : flag, =-1 if convergence is not obtained   *
  !*                     =1 if convergence is ok.             *
  !*        Gamma: starting value for the eigenvalue as input *
  !*               approximation of the eigenvalue with preci-*
  !*               sion dta in output.                        *
  !*          X1 : contains in output the associated eigen-   *
  !*               vector.                                    *
  !*                                                          *
  !************************************************************
  Subroutine IIM(eps, dta, m, n, A, it, gamma, X1)
    integer, parameter :: NMAX = 500
    real*8 A(NMAX,NMAX), X1(NMAX)
    real*8 eps, dta, gamma
    real*8  p0,phi,s,t0
    real*8 W(NMAX), X0(NMAX)
    integer LP(NMAX)

    do i=1, n
	  A(i,i) = A(i,i) - gamma
    end do
    do k=1, n-1
      p0=A(k,k); l0=k
      do i=k+1, n
		if (dabs(A(i,k)) > dabs(p0)) then
	      p0=A(i,k); l0=i
		end if
      end do 
      LP(k)=l0
      if (dabs(p0) < eps) then
	    p0=eps; A(l0,k)=eps
      end if
      if (l0.ne.k) then
	    do j=k, n
	      t0=A(k,j); A(k,j)=A(l0,j); A(l0,j)=t0
		end do
      end if
	  do i=k+1, n
	    A(i,k)=A(i,k)/p0
	    do j=k+1, n
	      A(i,j)=A(i,j)-A(i,k)*A(k,j)
        end do
      end do
	end do !k loop

    if (dabs(A(n,n)) < eps)  A(n,n)=eps
    do i=1, n
	  X0(i)=1.d0/sqrt(1.d0*i)
    end do 
    it=-1; l=1
    do while (it==-1.and.l<=m)
      do i=1, n
	    W(i)=X0(i)
      end do
      do k=1, n-1
        l0=LP(k)
        if (l0.ne.k) then
          t0=W(k); W(k)=W(l0); W(l0)=t0
        end if
        do i=k+1, n
		  W(i)=W(i)-A(i,k)*W(k)
        end do
      end do
      X1(n)=W(n)/A(n,n)
      do i=n-1, 1, -1
        s=0.d0
        do j=i+1, n
		  s = s + A(i,j)*X1(j)
        end do
	    X1(i)=(W(i)-s)/A(i,i)
      end do
      p0=0.d0
      do i=1, n
        if (dabs(X1(i)) > dabs(p0))  p0=X1(i)
      end do
      do i=1, n
        X1(i)=X1(i)/p0
      end do
      phi=0.d0
      do i=1, n
        s=dabs(X1(i)-X0(i))
        if (s > phi)  phi=s
      end do

      if (phi < dta) then
        gamma = gamma + 1.d0/p0
        it=1
      else
        do i=1, n
		  X0(i)=X1(i)
        end do 
        l = l + 1
      end if
	end do !while
	return
  End  !IIM  	  						    


  !******************************************************
  !* INPUTS:                                            *
  !* EPS : precision (Double)                           *
  !* D1  : precision d1 (Double)                        *
  !* D2  : precision d2 (Double)                        *
  !* M   : maximum number of iterations (integer)       *
  !* N   : order of matrix A (integer)                  *
  !* A   : input matrix to study (of MAT type)          *
  !* -------------------------------------------------- *
  !* OUTPUTS:                                           *
  !* IT  : -1 if no convergence is obtained (integer)   *
  !* R   : table of eigenvalues (of VEC type)           *
  !* VX  : table of eigenvectors (of MAT type)          *
  !******************************************************
  Subroutine EPMRI(eps, d1, d2, m, n, A, it, R, VX)
    integer, parameter :: NMAX = 500
    real*8 A(NMAX,NMAX), R(NMAX), VX(NMAX,NMAX)
    real*8 X(NMAX), A1(NMAX,NMAX)

    do i=1, n
      do j=1, n
        A1(i,j) = A(i,j)
      end do
    end do

    call VAMR(eps,d2,m,n,A1,it,R)

! restore A1 after VAMR
    do i=1, n
      do j=1, n
        A1(i,j) = A(i,j)
      end do
	end do   

    j=1
    do while (it==1.and.j<=n)

    call IIM(eps,d1,m,n,A1,it,R(j),X)
  
! restore A1 after IIM
      do i=1, n
        do k=1, n
          A1(i,k) = A(i,k)
        end do
        VX(i,j) = X(i)
      end do
      j = j + 1     
    end do !while
	
	return
  End

! End of file elpro.f90


!!if there is an external function, it will be here.
! Subroutine SelfDefined(matrix_1,matrix_2,rowNum_1,colNum_1,rowNum_2,colNum_2,matrix_result)


! end Subroutine SelfDefined

